function CellMaps(OriStat, labelimg, background, outdir, outformat, alpha, DI_threshold, showflag);

global forceGlobalAlpha;
global alpha;
global outformat;

if ~forceGlobalAlpha, alpha = 0.05; end

if nargin < 5
    outformat = '.jpg';
end

if nargin < 6
%     alpha = 0.05; % threshold for p_value

end

if nargin < 7
    DI_threshold = 0.33; % threshold for direction index.  
end

if nargin < 8
    showflag = 1;
end

stimStartShift = -45;  % in deg
oriPlotShift = -90;

% probability map is scaled between 10e-0 to 10e-2.

log_start = 0;
log_end = 2;

% scaling for tuning map 

narrowest = 15;
widest = 45;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ncells = size(OriStat);

red = [1, 0, 0];
blue = [0, 0, 1];

BlueRed_table = MakeColorGradationTable(blue, red, 100);

hue_table = MakeHueTable(100);

outdir_root = outdir;
outdir = [outdir, '\CellMap_'];
outformat2 = '';
global procroot

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% save background image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

imwrite(background, [outdir, 'avg', outformat]);

idx_resp = [OriStat.p_value_resp] <= alpha;
idx_sel = [OriStat.p_value_sel] <= alpha;

respIndices = find([OriStat.p_value_resp] > alpha); % indices of unresponsive cells
resplabelimg = labelimg;
for i = 1 : length([OriStat.p_value_resp])
    j = find( labelimg == i );
    if ~isempty( find( i == respIndices ) )
        resplabelimg(j) = 0;
    end
end
save( [outdir_root '_cells.mat' ], '-append', 'resplabelimg' );

selIndices = find([OriStat.p_value_sel] > alpha); % indices of unselective cells
sellabelimg = labelimg;
for i = 1 : length([OriStat.p_value_sel])
    j = find( labelimg == i );
    if ~isempty( find( i == selIndices ) )
        sellabelimg(j) = 0;
    end
end
save( [outdir_root '_cells.mat' ], '-append', 'sellabelimg' );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% probability maps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

logp_resp = (-log10([OriStat.p_value_resp]) - log_start) / (log_end - log_start);
logp_sel = (-log10([OriStat.p_value_sel]) - log_start) / (log_end - log_start);

MakeCellMap(logp_resp, 'all', labelimg, background, BlueRed_table, [outdir, 'p_value_resp', outformat], showflag);
MakeCellMap(logp_sel, 'all', labelimg, background, BlueRed_table, [outdir, 'p_value_sel', outformat], showflag);
MakePvalueLegend(alpha, log_start, log_end, procroot);



CellMapHistDensity([OriStat.p_value_resp], 'all', 0.005:0.01:0.995, [0 1], 'p', [outdir, 'p_value_resp_hist', outformat2]);
CellMapHistDensity([OriStat.p_value_sel], 'all', 0.005:0.01:0.995, [0 1], 'p', [outdir, 'p_value_sel_hist', outformat2]);

CellMapHistDensityDual([OriStat.p_value_sel], 'all', [OriStat.p_value_resp], 'all', 0.005:0.01:0.995, [0 1], 'p', [outdir, 'p_value_resp_sel_hist', outformat2]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% angle maps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ori_vector_angle = [OriStat.ori_vector_angle]/180;
MakeCellMap(ori_vector_angle, idx_resp, labelimg, background, hue_table, [outdir, 'ori_vector_angle', outformat], showflag);
CellMapHistDensity(-1 * (mod(ori_vector_angle * 180 - stimStartShift, 180) + oriPlotShift), idx_resp, -89:2:89, [-90 90], 'deg', [outdir, 'ori_vector_angle_hist', outformat2]);

dir_vector_angle = [OriStat.dir_vector_angle]/360;
dir_vector_angle_idx = ~([OriStat.p_value_sel] > alpha | [OriStat.DI] < DI_threshold);
MakeCellMap(dir_vector_angle, dir_vector_angle_idx, labelimg, background, hue_table, [outdir, 'dir_vector_angle', outformat], showflag);

best_dir_dp = [OriStat.best_dir_dp]/360;
best_dir_dp_idx = ~([OriStat.p_value_sel] > alpha | [OriStat.DI_dp] < DI_threshold);
MakeCellMap(best_dir_dp, best_dir_dp_idx, labelimg, background, hue_table, [outdir, 'best_dir_dp', outformat], showflag);

best_dir_interp = [OriStat.best_dir_interp]/360;
best_dir_interp_idx = ~([OriStat.p_value_sel] > alpha | [OriStat.DI_interp] < DI_threshold);
MakeCellMap(best_dir_interp, best_dir_interp_idx, labelimg, background, hue_table, [outdir, 'best_dir_interp', outformat], showflag);

best_dir_fit = [OriStat.best_dir_fit]/360;
best_dir_fit_idx = ~([OriStat.p_value_sel] > alpha | [OriStat.DI_fit] < DI_threshold);
MakeCellMap(best_dir_fit, best_dir_fit_idx, labelimg, background, hue_table, [outdir, 'best_dir_fit', outformat], showflag);

best_ori_fit = [OriStat.best_ori_fit]/180;
MakeCellMap(best_ori_fit, idx_resp, labelimg, background, hue_table, [outdir, 'best_ori_fit', outformat], showflag);
CellMapHistDensity(-1 * (mod(best_ori_fit * 180 - stimStartShift, 180) + oriPlotShift), idx_resp, -89:2:89, [-90 90], 'deg', [outdir, 'best_ori_fit_hist', outformat2]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   DI maps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DI = [OriStat.DI];
MakeCellMap(DI, idx_sel, labelimg, background, BlueRed_table, [outdir,'DI',outformat], showflag);

OI = [OriStat.OI];
MakeCellMap(OI, idx_resp, labelimg, background, BlueRed_table, [outdir,'OI',outformat], showflag);

DI_dp = [OriStat.DI_dp];
MakeCellMap(DI_dp, idx_sel, labelimg, background, BlueRed_table, [outdir,'DI_dp',outformat], showflag);

DI_interp = [OriStat.DI_interp];
MakeCellMap(DI_interp, idx_sel, labelimg, background, BlueRed_table, [outdir,'DI_interp',outformat], showflag);

DI_fit = [OriStat.DI_fit];
MakeCellMap(DI_fit, idx_sel, labelimg, background, BlueRed_table, [outdir,'DI_fit',outformat], showflag);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tuning width and reponse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dir_tuning_width = (widest - [OriStat.dir_tuning_width]) / (widest - narrowest); 
MakeCellMap(dir_tuning_width, idx_sel, labelimg, background, BlueRed_table, [outdir, 'dir_tuning_width', outformat], showflag);

ori_tuning_width = [OriStat.ori_tuning_width];
ori_tuning_width2 = (widest - ori_tuning_width) / (widest - narrowest); 
MakeCellMap(ori_tuning_width2, idx_resp, labelimg, background, BlueRed_table, [outdir, 'ori_tuning_width', outformat], showflag);
CellMapHistDensity(ori_tuning_width, idx_resp, 1:2:89, [0 90], 'deg', [outdir, 'ori_tuning_width_hist', outformat2]);

ori_vector_tune = [OriStat.ori_vector_tune];
CellMapHistDensity(ori_vector_tune, idx_resp, 0.005:0.01:0.995, [0 1], '|S|', [outdir, 'ori_vector_tune_hist', outformat2]);

R_best_ori = [OriStat.R_best_ori];
CellMapHistDensity(R_best_ori, idx_resp, 0.005:0.01:0.495, [0 0.5], 'ratio change', [outdir, 'R_best_ori_hist', outformat2]);

R_min_ori = [OriStat.R_min_ori];
CellMapHistDensity(R_min_ori, idx_resp, 0.005:0.01:0.495, [0 0.5], 'ratio change', [outdir, 'R_min_ori_hist', outformat2]);

% CellMapHistDensityDual(R_min_ori, idx_resp, R_best_ori, idx_resp, 0.005:0.01:0.495, [0 0.5], 'ratio change', [outdir, 'R_min_best_ori_hist', outformat2]);

ori_A1 = [OriStat.ori_A1];
CellMapHistDensity(ori_A1, idx_resp, 0.005:0.01:0.495, [0 0.5], 'ratio change', [outdir, 'ori_A1_hist', outformat2]);

ori_vector_angle = [OriStat.ori_vector_angle]/180;
ori_vector_angle([OriStat.p_value_resp] > alpha) = NaN;
R_best_ori = [OriStat.R_best_ori];
R_best_ori([OriStat.p_value_resp] > alpha) = NaN; % this will be the lightness
ori_vector_tune = [OriStat.ori_vector_tune];
ori_vector_tune([OriStat.p_value_resp] > alpha) = NaN; % this will be the saturation
hls.th = background;
hls.th(:) = 0;
hls.max_change = background;
hls.max_change(:) = 0;
hls.tune = background;
hls.tune(:) = 0;
for i = 1 : length(R_best_ori)
    k = find( labelimg == i );
    hls.th(k) = ori_vector_angle(i);
    hls.max_change(k) = R_best_ori(i);
    hls.tune(k) = ori_vector_tune(i);
end
WriteHlsMap( hls, [outdir,'ori_fit_'], max(max(hls.max_change)), max(max(hls.tune)) );
% WriteHlsMap( hls, [outdir,'ori_fit_'], max(max(hls.max_change)), 0.4 );
% WriteHlsMap( hls, [outdir,'ori_fit_'], 0.2, 0.4 );
max(max(hls.tune))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% scatter plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [0 1.1 * max(R_best_ori)] - for some reason this doesn't work for pooled
% OriStat files... going to just go for the obvious and make fixed XLim. - leongjcs 070129

CellMapScatter(R_min_ori, R_best_ori, idx_resp, 1, [0.05:0.1:0.45], [0 0.45], ...
    'min ratio change', [0 0.45], 'best ratio change', [outdir, 'R_min_best_scatter', outformat2]);

CellMapScatter(-1 * (mod(ori_vector_angle * 180 - stimStartShift, 180) + oriPlotShift), ori_vector_tune, ...
    idx_resp, 1, [-90 -45 0 45 90], [-90 90], 'interp best angle [deg]', ...
    [0 1.05], 'tuning width [ratio]', [outdir, 'scatter_vector_tune_vs_vector_angle', ...
    outformat2]);

close all;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CDF plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

binCenters = [0:180/4:180];
bins = mod(binCenters - 22.5, 180);
colors = binCenters / 180;
values = BinOriStat( OriStat, 'ori_vector_angle', bins, 'ori_vector_tune');
CellMapMultiCDFPlot( values, colors, hue_table, 'OSI', [outdir, 'CDF_ori_vector_tune'] );